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Abstract 

We develop novel methods to compute auto-correlation functions, or power spectral densi- 
ties, for chaotic dynamical systems generated by an inverse method whose starting point is an 
invariant distribution and a two-form. In general, the inverse method makes some aspects of 
chaotic dynamics calculable by methods familiar in quantum field theory. This approach has the 
numerical advantage of being amenable to Monte-Carlo parallel computation. We demonstrate 
the approach on a specific example, and show how auto-correlation functions can be computed 
without any direct numerical simulation, by Pade approximants of a short time expansion. 
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1 Introduction 



Due to their sensitive dependence on initial conditions, the precise behavior of a chaotic dynam- 
ical system over long times is unpredictable in practice, even though the underlying system is 
deterministic. However there is no fundamental obstruction to calculating statistical properties 
of the long time behavior of chaotic systems. 

Typically, statistical properties of chaotic systems are calculated by direct numerical simu- 
lation, with the assumption that the statistics, in the form of time averages, converge quickly 
compared to the rate at which numerical errors accumulate due to the sensitivity to initial con- 
ditions. Direct numerical simulation of systems with a very large number of degrees of freedom 
requires significant computational power which in many cases is not yet available. Moreover, it 
is difficult to gain theoretical insight from direct simulations. 

In a previous paper pQ, which we review in section II, we described an alternative approach 
to direct numerical simulation. This approach involves an inverse method whereby chaotic 
dynamical systems are generated, given a two-form and exact statistical information given by 
an invariant distribution over phase space. Several examples of chaotic systems were given for 
which there was precise agreement between equal time moments calculated by direct numerical 
simulation or by using the initial invariant distribution. 

There are several advantages to the inverse approach over direct simulation. In principle, the 
inverse approach can be applied to systems with a large number of degrees of freedom, without 
placing intense demands on computational power. It is also possible to begin a classification 
of chaotic systems in terms of analytic properties of the statistics. We expect that, using the 
inverse approach, many theoretical insights into properties of chaotic dynamical systems are 
forthcoming. 

In pQ, only equal time moments (and cumulants) were computed for chaotic systems gen- 
erated by the inverse approach. The present work is a sequel to pQ, in which we will apply 
the inverse method to the computation of auto-correlation functions, whose Fourier transform 
is the power spectral density. Ordinarily, the auto-correlation or power spectrum is computed 
by a direct numerical simulation with a long run-time. However, for chaotic systems generated 
by the inverse method, the auto-correlation can be computed by a Monte-Carlo simulation, 
involving many fixed time simulations from different initial conditions, as opposed to one long 
simulation. This has the great advantage of being amenable to parallel computation. We discuss 
the Monte-Carlo approach in section III. Although we only consider auto-correlation functions 
here, we emphasize that our approach is applicable to any correlation functions of phase space 
variables at different times. 

For chaotic systems generated by the inverse method, it is also possible to calculate auto- 
correlation functions without direct numerical simulation. In section IV, we show that the 
inverse approach yields a small t expansion for the auto-correlation G(t), whose validity can 
be extended deeper into the complex t plane by Pade approximants. We shall see that this 
approach shows good agreement with results obtained by direct numerical simulation. 

Power spectra of dynamical systems are sometimes obtained from Fourier decomposition of 
a signal within some time window. For a chaotic system, this approach is rife with ambiguities 
and pitfalls, unlike the Fourier transform of the auto-correlation function which is well defined. 
However these pitfalls are themselves interesting and can be well understood using the inverse 
approach, as we shall see in section V. 
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2 Review of the inverse approach 



We give a very brief review of the inverse method here. For details the reader is referred to [T]. 
We consider deterministic dynamical systems, with phase space coordinates x\(t) ■ ■ ■ xjsr(t) and 
equations of motion 

d/Xfi , , , . 

-— = v n {x 1 ,---x N ) (2.1) 

A distribution over initial conditions p[x] which is left invariant by the time evolution is a solution 
of the zero diffusion limit of the static Fokker-Planck equation, 

V • (pv) = , (2.2) 

For N = 3, this implies that pv can be written as the curl of a vector field, 

-. V x A , s 

v = , 2.3 

P 

and the inverse method amounts to choosing A and p, such that v is polynomial in the phase 
space coordinates x n . For N > 3, it is convenient to use the language of differential forms, in 
which case (12.11) becomes 



d*(pv) = (2.4) 

where v is the velocity one form v m dx m , d is the exterior derivative, and * indicates the Hodge 
dual. This implies that 

* dA (o 
v = (2.5) 

P 

where A is an N — 2 form. The problem is then to choose p and A such that v is polynomial in 
x n . This was done in pQ for a few increasingly complex analytic structures of p and A. 

Since *A is a two-form, which is more conveniently written down for large N than A, it can 
be said that an invariant distribution and a two-form uniquely specify dynamical systems. This 
is very similar to the specification of Hamiltonian dynamical systems by a Hamiltonian and a 
symplectic form, although far more general. For example N need not be even, and a Liouville 
theorem, d*v = 0, need not apply. There also are not necessarily any conserved quantities. 

Amongst the simplest distributions, so far as analytic structure is concerned, are polynomial 
distributions for which the real zeroes form a closed manifold inside of which the polynomial is 
positive. For these, 

A = p 2 n (2.6) 

were p is a polynomial, and f2 is a polynomial N — 2 form, so that 

v = *(pdQ + 2dp AQ) (2.7) 

Chaotic dynamical systems of this type are invariably repellers (see [T]). Chaotic attractors 
arise from distributions and two-forms with more complicated analytic structure, which are also 
discussed in |Tj. 
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Although a distribution p satisfying (12, lj) is an invariant distribution over initial conditions, 
it does not necessarily describe the statistics of a single chaotic trajectory. It must be verified 
that the dynamics is chaotic and that p, or strictly speaking its projection, is an ergodic measure, 
having no convex decomposition into independent invariant measures p ^ xp\ + (1 — x)pi. In 
many instances it is necessary to modify the domain of support of the initial distribution p to 
make it ergodic. For chaotic invariant sets with fractional information dimension, a function 
p(x), with support in N dimensions, can not be ergodic, although it may still contain exact 
information about the statistics, either via projection to lower integer dimension or, after Fourier 
transformation, as the generator of polynomial moments. The question of the general relevance 
of the distributions p arising in the inverse approach has not been definitively resolved. However, 
for the chaotic systems which have been considered so far in [1] , the distributions used to generate 
the dynamical systems give extremely accurate results for moments and cumulants (connected 
equal time correlation functions). Of course it is possible that these systems have information 
dimension which is very nearly N, in which case p yields, at worst, an extremely accurate 
approximation to the exact statistics. 

We emphasize that the results in this article are not restricted to chaotic dynamical systems 
generated by the inverse method. In fact, they may be applied to any system for which there is 
known information about the invariant measure over phase space or equal time moments. The 
point of this article is to use known information about static (time independent) statistics to 
facilitate the computation of time dependent quantities, such as auto-correlation functions and 
power spectral densities. 

3 Monte-Carlo computation of the auto-correlation 
and power spectrum 

Power spectra are useful signatures of dynamical systems derived from time series data. The 
power spectrum associated with a signal x(t) is defined as the squared amplitude of the Fourier 
transform of x(t). However, the Fourier transform of an a-periodic chaotic signal x(t) is not 
generally well defined since x(t) does not fall of as t — > ±oo. The power spectral density on the 
other hand is well defined; 

/oo 

dA exp(ic;A)G(A) (3.8) 
-oo 

where G(A) is the auto-correlation function, 
1 f T / 2 

G(A) = lim - / dtx(t)x(t + A) . (3.9) 

T-+oo T 

Due to the tendency of chaotic systems to 'forget' their initial conditions, G(A) falls off with 
large A (assuming < x >= 0) sufficiently rapidly that the Fourier transform (|3.8j) is well defined. 
The total power between frequencies ui and lo + dtu is given by G{oj)duj. 

A remarkable property of chaotic systems is that time averages of functions of phases space 
over a trajectory are equal to spatial averages with respect to an ergodic measure: 

</>= lim - f T dtf(x(t)) = 

T^oo 1 Jq 



dp(X)f(X) (3.10) 
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The auto-correlation of a phase space variable X{ can therefore be written as 

Gi(A) = J d»(X)f A (X) (3.11) 

where 

f A (X)=Xi(t = 0)xi(t = A), with Xi(0) = X { . (3.12) 

Thus, when the exact invariant measure is known, the auto-correlation and power spectral 
density can be computed by Monte-Carlo methods, which entail multiple fixed time simulations 
of the dynamical system from randomly generated initial conditions rather than a single very 
long simulation. This approach is amenable to efficient parallel computation schemes. 

4 An example 

Consider a specific chaotic system of the type (|2.7p (studied in PQ) with 

/0= l_ x 4_ y 2_ z 6 

n = (xyz)dy + {y 2 )dz . (4.13) 
yielding 

v x = 13z xy — 6y — xy + 2y + yx — 2yx — 2yz + y x , 

v y = 8xV , (4.14) 
v z = — 9x yz + yz — yz — y z . 

Initial conditions with p(x, y,z) > 0, y > and z > give rise to chaotic trajectories whose 
statistics are described by an invariant distribution p = p for p(x, y, z) > 0, y > 0, z > 0, up to 
a normalization, and p = outside this domain. 

Next we will evaluate the auto-correlation G X (A) =< x(t)x(t + A) > for these chaotic 
trajectories. The auto-correlation can be evaluated by direct numerical simulation, taking the 
time average of x(t)x(t + A) over a long run-time. Alternatively, it can be evaluated using 

^M, (EH; 

G X (A) = J dXdYdZ~p(X,Y,Z) f A (X,Y,Z) (4.15) 

where 

F A (X, Y, Z) = x(t = 0)x(t = A) (4.16) 

with 

x(t = 0) = X, y(t = 0) = Y, z{t = 0) = Z (4.17) 

Evaluating ()4.15|) by Monte Carlo simulation amounts to generating initial conditions X, Y, Z 
randomly according to the distribution p, and simulating the time evolution from each of these 
points over a duration A. The results of both the Monte-Carlo calculation and a direct long- 
duration numerical simulation of (|4.14p are shown in figure 1, showing extremely good agreement. 
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The example we have given is a small dimension system, with N = 3, so there is little 
computational advantage to the Monte-Carlo approach over a direct long duration simulation. 
The advantage comes at large N, for which the Monte-Carlo approach will be far faster, as it is 
amenable to parallel computation. 




Figure 1: Auto-correlation functions calculated by Monte-Carlo simulation (dashed line) 
and direct long-duration simulation (solid line). The dimensionless number j3 parameter- 
izes a variation of the dynamical system (14.141) . obtained from p = 1 — x A — y 2 — z 6 , Q = 

(xyz)dy + f3{y 2 )dz. 



5 Calculating the auto-correlation without direct 
numerical simulation. 

Using (|3.1ip . the auto-correlation function Gj(A) can be expanded as a power series in A. To 
do so let us first expand (|3,12p in A: 

F A (X) ~ a?i(0) (xi(0) + A±i(0) + ^A 2 xi(0) ■ ■ ■ J , with x(0) = X . (5.18) 

Note that Xi^^- can be written as a total time derivative for odd n: 

dxi d (I 2 \ d 3 Xi d I d 2 Xi 1 / dxA 2 ] 
Xi ^ = dt{2 Xi h Xi ^ = dt[ Xi l^-2{^t ' 6tC (5 - 19) 
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For even n, n = 2m, 



d x 4 d , x / \ Tn ( d Xi 



where we will not bother to specify g n , since total time derivatives vanish when averaged over a 
chaotic trajector}0. Thus 

oo 

G(A) = (F A ) ~ «- A2m ( 5 - 21 ) 



m=0 

1 / /d m x; 



2\ 



Using the equations of motion 



"A = Vl{x) 



the terms (d/dt) m Xi can be related to functions on phase space. For example 

d 2 Xi dvi \ dvi(x) 



so that 
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Since F&(X) can not have singularities on the real axis (the real solutions of the equations of 
motion are assumed to be non-singular), the series (|5.2ip has a finite radius of convergence. We 
propose to evaluate G(A) by Monte-Carlo integration to evaluate the coefficients a m , followed 
by a Pade resummation of (|5.21l) . 

For the dynamical system (14.141) , the moments a m are evaluated with respect to the measure 
dn(x) = dxdydzp(x) 

p = 1 - x 4 - y 2 - z 6 for x 4 + y 2 + z 6 < 1 , 

p = for x 4 + y 2 + z 6 > 1 , (5.23) 

We have calculated ctQ ■ ■ ■ as for the auto-correlation G X (A) using Monte-Carlo simulation, to- 
gether with a significant amount of algebra which was also carried out with a computer. The 
corresponding (4|4) Pade approximant is a ratio of polynomials of degree 4 whose first 9 Taylor- 
Mclaurin series coefficients are ao" ,,Q! 8- The (4|4) Pade approximant to G(A) is plotted in 
figure 2, along with the 9'th order Taylor Mclauren series result and the result of direct numeri- 
cal simulation by a long run of the dynamical system. Note that the 9-th order Taylor-Mclaurin 
series begins to differ markedly from the result of direct numerical simultation at A ~ 0.35, 
well before the first zero, whereas the (4|4) Pade approximant gives accurate result for much 
larger values of A, and is an acceptable approximation up to a neighborhood of the first zero, 
A ~ 2. Higher order Pade approximants are needed to obtain more real zeros, and to provide 



1 We are assuming a bounded system with no explicit time dependence in the equations of motion. 
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a good approximation of the exact result for larger values of A. This requires more computing 
power than we have presently applied to the problem. Nevertheless, the initial results are very 
encouraging, suggesting that the auto-correlation may be computed without a direct simulation 
of the dynamical system. 
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Figure 2: The auto-correlation calculated by long-run numerical simulation and by Monte- 
Carlo short run simulations (top), by 9'th order Monte-Carlo/Taylor-Mclauren series (bot- 
tom left) , and by Monte-Carlo/ (4 1 4) Pade approximant (bottom right). The Pade result 
is a smooth continuous curve, with the exception of a very small neighborhood of the 
point A = 0.6529, at which there is a spurious pole. 

An oft stated property of chaotic power spectra is that there is non-zero exponentially small 
component at high frequency [3j[2], G(oS) ~ exp(— aw). The time-scale a is determined by the 
proximity of the nearest singularity of the auto-correlation G(A) to the real A axis. Note that 
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there can not be any singularities on the real A axis, as it is assumed that the time evolution 
of the dynamical system does not encounter any singularities. Due to the tendency of chaotic 
systems to 'forget' their intial conditions, one expects singularities of G(A) near the real axis 
to occur for small values of Re(A). This suggest that the high frequency behavior of the power 
spectral density can be extracted from the small A behavior of the auto-correlation. 

It may be possible to use the Pade approximant to the auto-correlation to get an estimate of 
the parameter a. The poles of the Pade approximant do not necessarily correspond to the true 
analytic structure, and may in general be spurious. In fact the (4|4) Pade approximant we have 
computed here has two poles which are likely spurious, including one on the positive real axis 
which is definitely spurious. These poles are extremely close to zeros of the Pade approximant. 
There is one pole which is not near any zero, at A = 1.241 Verifying that this yields a good 
first approximation to a requires consideration of higher orders in the Pade sequence, (n\n) for 
n > 4, which we have yet to attempt. 

6 Remarks on the Fourier transform of a chaotic sig- 
nal 

The Fourier transform of a chaotic trajectory x(t), 

/oo 
dte iujt Xi(t). (6.24) 
-oo 

is not well defined. Nevertheless, it is quite common to take the Fourier (or discrete Fourier) 
transform of chaotic time series with some window of fixed duration, squaring the amplitude to 
give a kind of power spectrum. 

Let us consider a function on phase space fuj{X), defined by 

r-2-w/u) 

M x ) = £j dtx ^ t)e ~ tuJt (6 - 25) 

where 

x(t = 0) = X (6.26) 

and x(t) satisfies the equations of motion. Next consider a chaotic trajectory X(t), and define 
< fui > t° be the time averaged value of fuj(X(t)) on this trajectory. The average can be 
computed by summing over values of f w {X(t)) at arbitrary regular time intervals At, yielding a 
result which is independent of the interval. If one chooses the interval At = 2tt/uj, one arrives 
at the formal result 

1 r T 

< f u >= lim - / dte-^Xiit) (6.27) 

T^oo 1 J Q 

In terms of the invariant measure d/j,(X) on the phase space of a chaotic system, ergodicity 
implies 



< 



f u >= J d^(X)UX) (6.28) 
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Note that different initial conditions x(0) which are different points on the same chaotic trajec- 
tory lead to a different overall phase in (|6.27p . On the other hand, (|6.28j) . too which (|6.27p is 
formally equivalent, has no reference to an initial condition. This is consistent only if < f u >= 0. 

It is interesting to try to compute (|6.28p by Monte-Carlo simulation. Note that any non- 
zero result is pure error, which we shall see is closely related to the power spectral density. The 
Monte-Carlo calculation amounts to generating a set of n initial conditions randomly, according 
to the distribution of the chaotic invariant set, and then running the dynamical system from 
each initial condition for a duration 2it/uj. Given the chaotic 'loss of information' about initial 
conditions with time, one expects this to yield a very similar result to direct numerical simulation 
from a single initial condition over a duration 2im/u. 

Calculating < f u > by Monte-Carlo simulation, any non-zero result is pure error, the size 
of which is related to the variance of f u . The variance of f u can be expressed in terms of the 
auto-correlation, 



< |/ w - < fu, > I 2 > 

x(t = 0)=X 



2tt/lu p2tt/uj 

dt / dt' I dfJL(X)x(t)x(t r )e iu3 ^ 
o Jo 



so that 



< \foj~ < > | 2 >= 



2n/u> /-2tt/uj 

dt / dt'Gtt-t'y^-^ 

ii Jo 



(6.29) 
(6.30) 

(6.31) 



Note that (|6.3ip is very similar to the power spectral density, with the exception of the boundaries 
of integration. 

We have calculated a 'power spectrum' for the dynamical system defined by (|4.13p . defined 
as the square of the amplitude of a discrete Fourier transform within some time window, or of 
a Monte-Carlo approximation to f u . The results of both calculations for (3 = 1 are plotted in 
figure 3. 





Figure 3: 'Power spectrum' computed by Fast Fourier Transform of a chaotic signal (on 
the left) and by Monte-Carlo simulation (on the right). 



Due to the tendency of chaotic systems to 'forget' their initial conditions, on might expect the 
Fourier transform of a run of finite but long duration to be qualitatively similar to the Monte- 
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Carlo result, which involves shorter runs from randomly generated initial conditions. Indeed 
the results plotted in figure 3 have some crude qualitative equivalence. In some sense, both 
these results can be viewed as nothing but finite total simulation time errors, the exact result 
vanishing identically, where the error is related to a quantity (|6.3ip similar, but not equivalent, 
to the power spectral density. 

7 Conclusions 

Direct numerical simulation of chaotic dynamical systems is a viable method to compute their 
statistics. However this approach rarely yields theoretical insight, and places extreme or pro- 
hibitive demands on computational resources for systems with a very large number of degrees 
of freedom. The intent of this work has been to develop an inverse method, introduced in [1], 
which permits the calculation of statistical quantities of chaotic systems without the use of direct 
numerical simulation over long time durations. 

The inverse method generates chaotic dynamical systems, given an invariant measure and a 
two-form. At present we we have tested the approach on systems with a small number of degrees 
of freedom, by computing static statistical quantities in [1] such as equal time correlations, and 
time dependent statistical quantities such as auto-correlation functions in the present article. 

Already at the level of a few degrees of freedom, one should be able to make theoretical 
progress which is not possible with studies involving only direct numerical simulation. In par- 
ticular it should be possible to construct a classification of chaotic dynamical systems according 
to the phase space analytic structure of an invariant distribution and a two-form. 

Given the invariant measure of a chaotic dynamical system, which is the starting point of the 
inverse approach, one can write closed form expressions for equal time correlation functions, for 
which analytic approximation schemes may exist, such as saddle point expansions. Non-equal 
time correlation functions, such as auto-correlations, may also be computable without recourse 
to direct numerical simulation, by Pade re-summation of short time expansions. In fact, many 
of the calculations arising using the inverse approach bear a close resemblance to calculations in 
quantum field theory, which may yield useful ideas in the present context of chaotic dynamics. 

The true power of the inverse method will come for a large number of degrees of freedom, since 
this approach is amenable to parallel computation. There are numerous interesting problems 
which remain to be solved, in particular that of reverse engineering a chaotic dynamical system 
of particular physical interest, such as a turbulent fluid. Efforts along these lines are underway. 

There remain subtle and interesting theoretical problems concerning the meaning of the 
invariant distribution, which is the starting point of the inverse method, when the information 
dimension is fractional. These subtleties have been discussed in pQ and remain to be resolved. 
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